home *** CD-ROM | disk | FTP | other *** search
/ C/C++ Users Group Library 1996 July / C-C++ Users Group Library July 1996.iso / listings / v_07_05 / v7n5110a.txt < prev   
Text File  |  1989-04-15  |  9KB  |  305 lines
  1. Listing 1
  2.  
  3.  
  4. /*  FFT - Cooley-Tukey Fast Fourier Transform
  5.  *
  6.  *  Purpose:   To calculate the discrete Fast
  7.  *             Fourier Transform (or its inverse)
  8.  *             of a set of complex data elements.
  9.  *
  10.  *  Setup:     COMPLEX data[];
  11.  *             int n;
  12.  *             int flag;
  13.  *             void fft();
  14.  *
  15.  *  Call:      fft(data, n, flag);
  16.  *
  17.  *  Where:     data is an array of n+1 data
  18.  *               elements, where data[0] is
  19.  *               unused.
  20.  *             n is the number of valid data
  21.  *               elements in "data" ( i.e. -
  22.  *               data[1...n] ). It MUST be
  23.  *               an integer power of two.
  24.  *             flag is a Boolean flag which
  25.  *               if zero indicates that the
  26.  *               inverse Fourier transform is
  27.  *               to be calculated; otherwise
  28.  *               the Fourier transform is
  29.  *               calculated.
  30.  *
  31.  *  Result:    The Fourier transform (or the
  32.  *             inverse Fourier transform) is
  33.  *             calculated in place and returned
  34.  *             in "data"; the original data
  35.  *             elements are overwritten.
  36.  *
  37.  *  Note:      The data type COMPLEX is:
  38.  *
  39.  *                  typedef struct complex
  40.  *                  {
  41.  *                    float real;
  42.  *                    float imag;
  43.  *                  }
  44.  *                  COMPLEX;
  45.  */
  46.  
  47. void fft(data, n, flag)
  48. COMPLEX data[];
  49. int n;
  50. int flag;
  51. {
  52.   int a, b, c, d, e; /* Counter variables */
  53.   COMPLEX temp;      /* Swap variable */
  54.   COMPLEX u, v;      /* Temporary variables */
  55.   double theta;
  56.   double PI = 3.14159265358979;
  57.  
  58.   /* Perform bit reversal */
  59.  
  60.   theta = flag ? PI : -PI;
  61.   b = n/2 + 1;
  62.   for (a = 2; a < n; a++)
  63.   {
  64.     if (a < b)
  65.     {
  66.       temp.real = data[a].real;
  67.       temp.imag = data[a].imag;
  68.       data[a].real = data[b].real;
  69.       data[a].imag = data[b].imag;
  70.       data[b].real = temp.real;
  71.       data[b].imag = temp.imag;
  72.     }
  73.     c = n/2;
  74.     while (c < b)
  75.     {
  76.       b -= c;
  77.       c /= 2;
  78.     }
  79.     b += c;
  80.   }
  81.  
  82.   /* Calculate Fourier transform */
  83.  
  84.   c = 1;
  85.   while (c < n)
  86.   {
  87.     d = c;
  88.     c *= 2;
  89.     u.real = 1.0;
  90.     u.imag = 0.0;
  91.     v.real = cos(theta/d);
  92.     v.imag = sin(theta/d);
  93.     for (b = 1; b <= d; b++)
  94.     {
  95.       for (a = b; a <= n; a += c)
  96.       {
  97.         e = a + d;
  98.         temp.real = data[e].real * u.real -
  99.                     data[e].imag * u.imag;
  100.         temp.imag = data[e].real * u.imag -
  101.                     data[e].imag * u.real;
  102.         data[e].real = data[a].real - temp.real;
  103.         data[e].imag = data[a].imag - temp.imag;
  104.         data[a].real += temp.real;
  105.         data[a].imag += temp.imag;
  106.       }
  107.       temp.real = u.real;
  108.       u.real = u.real * v.real -
  109.                u.imag * v.imag;
  110.       u.imag = temp.real * w.imag +
  111.                u.imag * v.real;
  112.     }
  113.   }
  114.  
  115.   /* Normalize if not inverse transform */
  116.  
  117.   if (flag)
  118.     for (a = 1; a <= n; a++)
  119.     {
  120.       data[a].real /= n;
  121.       data[a].imag /= n;
  122.     }
  123. }
  124.  
  125.  
  126. Listing 1 - Fast Fourier Transform (Cooley-Tukey)
  127. -------------------------------------------------
  128.  
  129.  
  130.  
  131.  
  132.  
  133. Listing 2
  134.  
  135.  
  136. /*  FWT - Shank's Fast Walsh Transform
  137.  *
  138.  *  Purpose:   To calculate the discrete Fast
  139.  *             Walsh Transform (or its inverse)
  140.  *             of a set of real data elements.
  141.  *
  142.  *  Setup:     float data[];
  143.  *             int n;
  144.  *             int flag;
  145.  *             void fwt();
  146.  *
  147.  *  Call:      fwt(data, n, flag);
  148.  *
  149.  *  Where:     data is an array of n+1 data
  150.  *               elements, where data[0] is
  151.  *               unused.
  152.  *             n is the number of valid data
  153.  *               elements in "data" ( i.e. -
  154.  *               data[1...n] ). It MUST be
  155.  *               an integer power of two.
  156.  *             flag is a Boolean flag which
  157.  *               if zero indicates that the
  158.  *               inverse Walsh transform is
  159.  *               to be calculated; otherwise
  160.  *               the Walsh transform is
  161.  *               calculated.
  162.  *
  163.  *  Result:    The Walsh transform (or the
  164.  *             inverse Walsh transform) is
  165.  *             calculated in place and returned
  166.  *             in "data"; the original data
  167.  *             elements are overwritten.
  168.  *
  169.  *  Note:      The Walsh function coefficients
  170.  *             returned in data are arranged in
  171.  *             increasing Paley (dyadic) order.
  172.  */
  173.  
  174. void fwt(data, n, flag)
  175. float data[];
  176. int n;
  177. int flag;
  178. {
  179.   int a, b, c, d, e; /* Counter variables */
  180.   double temp;       /* Swap variable */
  181.  
  182.   /* Perform bit reversal */
  183.  
  184.   b = n/2 + 1;
  185.   for (a = 2; a < n; a++)
  186.   {
  187.     if (a < b)
  188.     {
  189.       temp = data[a];
  190.       data[a] = data[b];
  191.       data[b] = temp;
  192.     }
  193.     c = n/2;
  194.     while (c < b)
  195.     {
  196.       b -= c;
  197.       c /= 2;
  198.     }
  199.     b += c;
  200.   }
  201.  
  202.   /* Calculate Walsh transform */
  203.  
  204.   c = 1;
  205.   while (c < n)
  206.   {
  207.     d = c;
  208.     c *= 2;
  209.     for (b = 1; b <= d; b++)
  210.       for (a = b; a <= n; a += c)
  211.       {
  212.         e = a + d;
  213.         temp = data[e];
  214.         data[e] = data[a] - temp;
  215.         data[a] += temp;
  216.       }
  217.   }
  218.  
  219.   /* Normalize if not inverse transform */
  220.  
  221.   if (flag)
  222.     for (a = 1; a <= n; a++)
  223.       data[a] /= n;
  224. }
  225.  
  226.  
  227. Listing 2 - Fast Walsh Transform (Shank)
  228. ----------------------------------------
  229.  
  230.  
  231.  
  232.  
  233. Figure 1
  234.  
  235.  
  236.                             |<-------- ONE PERIOD -------->|
  237.  
  238.                          +1 |-------------------------------
  239. WAL(0), PAL(0), CAL(0)    0 +-------------------------------->
  240.                          -1 |
  241.  
  242.                          +1 |--------------+
  243. WAL(1), PAL(1), SAL(1)    0 +--------------|----------------->
  244.                          -1 |              +----------------
  245.  
  246.                          +1 |-------+               +-------
  247. WAL(2), PAL(3), CAL(1)    0 +-------|---------------|-------->
  248.                          -1 |       +---------------+
  249.  
  250.                          +1 |-------+       +-------+      
  251. WAL(3), PAL(2), SAL(2)    0 +-------|-------|-------|-------->
  252.                          -1 |       +-------+       +-------
  253.  
  254.                          +1 |---+       +-------+       +---
  255. WAL(4), PAL(6), CAL(2)    0 +---|-------|-------|-------|---->
  256.                          -1 |   +-------+       +-------+
  257.  
  258.                          +1 |---+       +---+   +-------+
  259. WAL(5), PAL(7), SAL(3)    0 +---|-------|---|---|-------|---->
  260.                          -1 |   +-------+   +---+       +---
  261.  
  262.                          +1 |---+   +---+       +---+   +---
  263. WAL(6), PAL(5), CAL(3)    0 +---|---|---|-------|---|---|---->
  264.                          -1 |   +---+   +-------+   +---+
  265.  
  266.                          +1 |---+   +---+   +---+   +---+
  267. WAL(7), PAL(4), SAL(4)    0 +---|---|---|---|---|---|---|---->
  268.                          -1 |   +---+   +---+   +---+   +---
  269.  
  270.                          +1 |-+   +---+   +---+   +---+   +-
  271. WAL(8), PAL(12), CAL(4)   0 +-|---|---|---|---|---|---|---|-->
  272.                          -1 | +---+   +---+   +---+   +---+
  273.  
  274.                          +1 |-+   +---+   +-+ +---+   +---+
  275. WAL(9), PAL(13), SAL(5)   0 +-|---|---|---|-|-|---|---|---|-->
  276.                          -1 | +---+   +---+ +-+   +---+   +-
  277.  
  278.                          +1 |-+   +-+ +---+   +---+ +-+   +-
  279. WAL(10), PAL(15), CAL(5)  0 +-|---|-|-|---|---|---|-|-|---|-->
  280.                          -1 | +---+ +-+   +---+   +-+ +---+
  281.  
  282.                          +1 |-+   +-+ +---+ +-+   +-+ +---+
  283. WAL(11), PAL(14), SAL(6)  0 +-|---|-|-|---|-|-|---|-|-|---|-->
  284.                          -1 | +---+ +-+   +-+ +---+ +-+   +-
  285.  
  286.                          +1 |-+ +-+   +-+ +---+ +-+   +-+ +-
  287. WAL(12), PAL(10), CAL(6)  0 +-|-|-|---|-|-|---|-|-|---|-|-|-->
  288.                          -1 | +-+ +---+ +-+   +-+ +---+ +-+
  289.  
  290.                          +1 |-+ +-+   +-+ +-+ +-+ +---+ +-+
  291. WAL(13), PAL(11), SAL(7)  0 +-|-|-|---|-|-|-|-|-|-|---|-|-|-->
  292.                          -1 | +-+ +---+ +-+ +-+ +-+   +-+ +-
  293.  
  294.                          +1 |-+ +-+ +-+ +-+   +-+ +-+ +-+ +-
  295. WAL(14), PAL(9), CAL(7)   0 +-|-|-|-|-|-|-|---|-|-